#### This script will extract data from BRT daily Tmax preds for all WT sites for all days within a single year
#### Data will be bound into a single data frame, and quantile will be run to determine the 90th, 95th, and 99th percentiles
#### Both the summary data frame and percentiles data frame will be saved

### Necessary library

library('SDMTools')

#### Location of BRT tmax ASCII files 

tmax.dir = '/home/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/brtpredsfinal/maxgzip/'

### Need to create base.pos with all positions to intersect with their respective lat/longs
### First read in the DEM ASCII as the template

base.asc = read.asc('/home/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/STATIC/dem_WTplusbuffer_LatLong_WGS1984_250mres.asc')

### Identify all positions where data exists (n should equal 499085)

base.pos = as.data.frame(which(is.finite(base.asc), arr.ind=T))

### Convert all remaining finite data in base.asc to NA

base.asc[,] = NA

### Get lat longs for data extraction

base.pos$lat = getXYcoords(base.asc)$x[base.pos$row]
base.pos$long = getXYcoords(base.asc)$y[base.pos$col]

### Save base.pos as an object to be loaded in the run script

save(base.pos, base.asc, file=paste(tmax.dir,'baseascdata.Rdata',sep=''))

### Establish a PBS directory to write .sh files to

tmp.pbs = '/home/99/jc152199/tmp.pbs9/'
setwd(tmp.pbs)

### Location of run script

script.file = paste(tmax.dir,'BRT_Extract_By_Year_2run.r',sep='')

#### Need to batch each year out to a different node

for (yearx in c(1951:2011))

	{
	
	tname = yearx #Get the just the name of the ASCII file from the file list
	
	tname = paste(yearx,'.sh',sep='') #Use tname to create the name for the .sh file
	
	arg.yearx = paste('yearx="',yearx,'" ',sep='') # Define the argument for tfile 
	
	zz = file(tname,'w') # This command creates a file with the name 'tname' The 'w' command opens the file and writes the following lines of text into it
		
		cat('##################################\n',file=zz) # First blank line of .sh file ### Note the use of '\n', any text following this will be on the next line
		cat('#!/bin/sh\n',file=zz) # First command line for HPC
		cat('cd $PBS_O_WORKDIR\n',file=zz) # Sets the directory to write out summary files at the above directory 'tmp.pbs'
		cat("R CMD BATCH --no-save --no-restore '--args ",arg.yearx,"' ",script.file,' ',gsub('.sh','',tname),'.Rout \n',sep='',file=zz) # This line tells R not to save to the .Rdata file and also not to restore the environment from the .Rdata file, very important when batching or you will get multiple fuckups, also defines the arguments, says where script to run is located, and defines the name for the R.out file
		cat('##################################\n',file=zz) # Last blank line of .sh file
	
	close(zz) # Closes the .sh file that you've just written
	
	# The follow line submits the job to the HPC, leave it commented out until you're ready to RUN all the jobs\
	# Always run the script first without removing the comment, this will still generate the .sh files for you to check before running
	# Once you're happy with that then submit a single .sh file and make sure it runs through to completion before batching out the rest with the following line
		
	system(paste('qsub -m n -l nodes=1:ppn=8 ',tname,sep=''))
	
	system(paste('sleep 1'))
	
	}
